\documentclass{article}
\usepackage{lmodern}
\usepackage[a4paper,scale=0.8]{geometry}
\usepackage{listings}
\usepackage[dvipsnames]{xcolor}
\lstset{
    language=Python,
    basicstyle=\ttfamily,
    breaklines=true,
    keywordstyle=\bfseries\color{NavyBlue}, 
    commentstyle=\itshape\color{black!50!white},
    stringstyle=\bfseries\color{PineGreen!90!black},
    columns=flexible,
    numbers=left,
    numbersep=2em,
    numberstyle=\footnotesize,
    frame=single,
    framesep=1em
}

\usepackage{pgf-umlcd}

\usepackage{ctex}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{float}
\usepackage{datetime}


\title{数值分析作业5}
\author{陈乐瑶  3210103710}
\date{\today}

\begin{document}
\maketitle

\part{A and B}
\begin{itemize}
    \item 线性方程组求解：np.linalg.solve ()
    \item 矩阵条件数计算：np.linalg.cond ()
\end{itemize}

\section{DLS via normal equations}
Normal Equation: 
\begin{align}
    G\alpha=c
\end{align}

实现代码如下：
\begin{lstlisting}
def DLS_via_normal_equation(x, y):
    if np.size(x) != np.size(y):
        print('ERROR: the size of two vectors does not match!')
        return
    G = np.zeros((3, 3))
    c = np.zeros(3)
    for i in range(3):
        for j in range(3):
            k = i+j
            G[i][j] = sum(np.power(x, k))
    for i in range(3):
        c[i] = sum(y*np.power(x, i))
    a = np.linalg.solve(G, c)
    print('cond(G): ', np.linalg.cond(G))
    return a
\end{lstlisting}

\section{DLS via QR factorization}
QR factorization:
\begin{align}
    \min \limits_{x \in \mathbb{R} ^n}\|Ax-b\|_2^2 &\Leftrightarrow  R_1 x = c,\\
    A &= QR,\\
    Q^T A = R = \begin{bmatrix}
        R \\
        O \\
    \end{bmatrix} &, \quad 
    Q^T b = \begin{bmatrix}
        c \\
        r \\
    \end{bmatrix}
\end{align}

实现代码如下：
\begin{lstlisting}
def DLS_via_QR_factorization(x, y):
    if np.size(x) != np.size(y):
        print('ERROR: the sizes of two vectors do not match!')
        return
    n = np.size(x)
    A = np.zeros((n, 3))
    b = np.zeros(n)
    for i in range(n):
        for j in range(3):
            A[i][j] = pow(x[i], j)
        b[i] = y[i]
    Q, R = np.linalg.qr(A)
    c = np.dot(np.transpose(Q), b)
    c = c[0:3]
    a = np.linalg.solve(R[0:3], c)
    print('cond(R): ', np.linalg.cond(R))
    return a
\end{lstlisting}

三组数据的二次多项式拟合结果如下：
\begin{figure}[htbp]
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5cm]{./pic/DLS_for_data_1.png}
        \caption{data1}
    \end{minipage}
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5cm]{./pic/DLS_for_data_2.png}
        \caption{data2}
    \end{minipage}
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5cm]{./pic/DLS_for_data_3.png}
        \caption{data3}
    \end{minipage}
    \caption{DLS via normal equations and QR factorization}
\end{figure}

矩阵条件数和多项式系数输出如下：
\begin{figure}[htbp]
    \includegraphics[width=14cm]{./pic/problemAB_output.png}
    \caption{output}
\end{figure}

\part{C}
\section*{(1)}
由几何级数 $\sum_{i=1}^{j} (-1)^{i-1} x^{i-1} = \frac{1-(-1)^j x^j}{1+x} $并对两边同乘$ (-1) ^j$可推知，$\frac{x^j}{1+x} = \sum_{i=1}^{j} (-1)^{i+j} x^{i-1} - \frac{(-1)^{j+1}}{1+x}$。
当$f(x) = \frac{1}{1+x}$时， (5.36) 的右式可写为$\int_{0}^{1} \frac{x^j}{1+x} \,dx = \int_{0}^{1} ( \sum_{i=1}^{j} (-1)^{i+j} x^{i-1} - \frac{(-1)^{j+1}}{1+x}) \,dx $。\\
由于$(-1)^{i+j}x^{i-1}$连续且$\sum_{i=1}^{j} (-1)^{i+j}x^{i-1}$一致收敛，
所以$\int_{0}^{1} \sum_{i=1}^{j} (-1)^{i+j} x^{i-1} \,dx = \sum_{i=1}^{j} (-1)^{i+j} \int_{0}^{1} x^{i-1} \,dx = \sum_{i=1}^{j} (-1)^{i+j} \frac{1}{i} , i \neq 0$。\\
而另一方面，$\int_{0}^{1} \frac{(-1)^{j+1}}{1+x} \,dx = (-1)^{j+1} \ln2 $。因此可推出 (5.36) $r_j = (-1)^j ln2 + \sum_{i=1}^{j} (-1)^{i+j} \frac{1}{i}$，其中易知，$(-1)^j$和$\sum_{i=1}^{j} (-1)^{i+j}$均为有理数，我们将这两个有理数记作$r_1$和$r_2$。\\
由 (5.35) 可知，$r = r_1 \ln2 + r_2$ 为 $\alpha$的线性组合，所以我们可以将$\alpha$表达为两个部分$\beta$和$\gamma$，即$\alpha = \beta \ln2 + \gamma$。
则有$H \beta = r_1, H \gamma = r_2$，下证$\beta$和$\gamma$均为有理数。\\
由于Hilbert矩阵正定可逆，且其行列式和伴随矩阵的元素均为原矩阵元素的四则运算结果，所以$H$的逆矩阵$H^{-1}$存在且其元素也均为有理数。
那么$\beta = H^{-1}r_1, \gamma = H^{-1}r_2$的各项均为有理数的四则运算，仍为有理数。另一方面，我们也可以通过高斯消元证明其为有理数的线性运算结果，仍为有理数。\\
从而我们可以得到 (5.35)解的表达式：$\forall j = 0, 1, \ldots, n, \alpha_j = \beta_j \ln2 + \gamma_j$，其中$\beta_j$和$\gamma_j$为有理数。

\section*{(2)}
用分数进行高斯消去，代码如下：
\begin{lstlisting}
def gauss_elimination(matrix, vector):
    if len(matrix) != len(vector):
        print('ERROR: the size of matrix and vector does not match!')
        return
    M = np.column_stack([matrix, vector])
    n = len(matrix)
    for i in range(n-1):
        for j in range(i+1, n):
            factor = Fraction(M[j, i], M[i, i])
            M[j, i:] -= factor * M[i, i:]
    x = np.zeros(n, Fraction)
    for i in range(n-1, -1, -1):
        x[i] = (M[i, n] - np.dot(M[i, i:n], x[i:n])) / M[i, i]
    return x
\end{lstlisting}

将（5.36）$r_j$的表达式分为两部分，$r_{j_1} = (-1)^j$以及$r_{j_2} = \sum_{i=1} ^j (-1)^{i+j} \frac{1}{j}$，从而$r_j = r_{j_1} \ln2 + r_{j_2}$。
那么只需要分别对增广矩阵$\left[H|r_{j_1}\right]$和$\left[H|r_{j_2}\right]$进行高斯消去即可得到$\beta_j$和$\gamma_j$.

\section*{(3) (4) (5)}
\begin{itemize}
    \item 用公式计算：$\alpha_j = \beta_j \ln2 + \gamma_j$，其中$\ln2$用机器精度逼近。
    \item 线性方程组计算：$H\alpha = r$，其中$r_j$用浮点数计算，$\ln2$用0.69315表示。
    \item 线性方程组计算：$(kI + H^*H)\alpha = H*r$，其中k为Tikhonov正则化参数，$H^*$为$H$的共轭转置。
\end{itemize}


不同n下的$\beta$，$\gamma$和Hilbert矩阵条件数的计算结果如下：

\begin{figure}[htbp]
    \includegraphics[width=16cm]{./pic/problemC_output_1.png}
    \caption{cond (H), $\beta$ and $\gamma$}
\end{figure}

三种情况计算得到的系数向量$\alpha$如下：

\begin{figure}[htbp]
    \includegraphics[width=16cm]{./pic/problemC_output_2.png}
    \caption{$\alpha$}
\end{figure}

以$\alpha_0, \alpha_1, \alpha_2$为例作图显示$\alpha$相对于泰勒展开系数的误差收敛性：

\begin{figure}[htbp]
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5cm]{./pic/convergence_of_a0.png}
        \caption{$\alpha_0$}
    \end{minipage}
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5cm]{./pic/convergence_of_a1.png}
        \caption{$\alpha_1$}
    \end{minipage}
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5cm]{./pic/convergence_of_a2.png}
        \caption{$\alpha_2$}
    \end{minipage}
    \caption{convergence of $\alpha_0, \alpha_1, \alpha_2$}
\end{figure}

三种$\alpha$的收敛性解释：
\begin{itemize}
    \item 用公式表达式计算时由于使用分数表示有理数进行计算，所以几乎不存在舍入误差。
    \item 由于Hilbert矩阵的条件数随着n的增大而指数增大，所以在n较大时，用浮点数进行矩阵计算时微小的舍入误差也会导致结果$\alpha$的大幅度误差。
    \item 而用Tikhonov正则化参数k来控制矩阵条件数，可以使得问题的解变得更加平滑，从而减小由于病态问题引起的误差。
\end{itemize}

\end{document}